home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qags.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  12.9 KB  |  571 lines
  1. /* integration/qags.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_errno.h>
  23. #include <gsl/gsl_integration.h>
  24.  
  25. #include "initialise.c"
  26. #include "set_initial.c"
  27. #include "qpsrt.c"
  28. #include "util.c"
  29. #include "reset.c"
  30. #include "qpsrt2.c"
  31. #include "qelg.c"
  32. #include "positivity.c"
  33.  
  34. static int qags (const gsl_function * f, const double a, const double
  35.   b, const double epsabs, const double epsrel, const size_t limit,
  36.   gsl_integration_workspace * workspace, double *result, double *abserr,
  37.   gsl_integration_rule * q);
  38.  
  39. int
  40. gsl_integration_qags (const gsl_function *f,
  41.               double a, double b,
  42.               double epsabs, double epsrel, size_t limit,
  43.               gsl_integration_workspace * workspace,
  44.               double * result, double * abserr)
  45. {
  46.   int status = qags (f, a, b, epsabs, epsrel, limit,
  47.                      workspace, 
  48.                      result, abserr, 
  49.                      &gsl_integration_qk21) ;
  50.   return status ;
  51. }
  52.  
  53. /* QAGI: evaluate an integral over an infinite range using the
  54.    transformation
  55.  
  56.    integrate(f(x),-Inf,Inf) = integrate((f((1-t)/t) + f(-(1-t)/t))/t^2,0,1)
  57.  
  58.    */
  59.  
  60. static double i_transform (double t, void *params);
  61.  
  62. int
  63. gsl_integration_qagi (gsl_function * f,
  64.               double epsabs, double epsrel, size_t limit,
  65.               gsl_integration_workspace * workspace,
  66.               double *result, double *abserr)
  67. {
  68.   int status;
  69.  
  70.   gsl_function f_transform;
  71.  
  72.   f_transform.function = &i_transform;
  73.   f_transform.params = f;
  74.  
  75.   status = qags (&f_transform, 0.0, 1.0, 
  76.                  epsabs, epsrel, limit,
  77.                  workspace,
  78.                  result, abserr,
  79.                  &gsl_integration_qk15);
  80.  
  81.   return status;
  82. }
  83.  
  84. static double 
  85. i_transform (double t, void *params)
  86. {
  87.   gsl_function *f = (gsl_function *) params;
  88.   double x = (1 - t) / t;
  89.   double y = GSL_FN_EVAL (f, x) + GSL_FN_EVAL (f, -x);
  90.   return (y / t) / t;
  91. }
  92.  
  93.  
  94. /* QAGIL: Evaluate an integral over an infinite range using the
  95.    transformation,
  96.    
  97.    integrate(f(x),-Inf,b) = integrate(f(b-(1-t)/t)/t^2,0,1)
  98.  
  99.    */
  100.  
  101. struct il_params { double b ; gsl_function * f ; } ;
  102.  
  103. static double il_transform (double t, void *params);
  104.  
  105. int
  106. gsl_integration_qagil (gsl_function * f,
  107.                double b,
  108.                double epsabs, double epsrel, size_t limit,
  109.                gsl_integration_workspace * workspace,
  110.                double *result, double *abserr)
  111. {
  112.   int status;
  113.  
  114.   gsl_function f_transform;
  115.   struct il_params transform_params  ;
  116.  
  117.   transform_params.b = b ;
  118.   transform_params.f = f ;
  119.  
  120.   f_transform.function = &il_transform;
  121.   f_transform.params = &transform_params;
  122.  
  123.   status = qags (&f_transform, 0.0, 1.0, 
  124.                  epsabs, epsrel, limit,
  125.                  workspace,
  126.                  result, abserr,
  127.                  &gsl_integration_qk15);
  128.  
  129.   return status;
  130. }
  131.  
  132. static double 
  133. il_transform (double t, void *params)
  134. {
  135.   struct il_params *p = (struct il_params *) params;
  136.   double b = p->b;
  137.   gsl_function * f = p->f;
  138.   double x = b - (1 - t) / t;
  139.   double y = GSL_FN_EVAL (f, x);
  140.   return (y / t) / t;
  141. }
  142.  
  143. /* QAGIU: Evaluate an integral over an infinite range using the
  144.    transformation
  145.  
  146.    integrate(f(x),a,Inf) = integrate(f(a+(1-t)/t)/t^2,0,1)
  147.  
  148.    */
  149.  
  150. struct iu_params { double a ; gsl_function * f ; } ;
  151.  
  152. static double iu_transform (double t, void *params);
  153.  
  154. int
  155. gsl_integration_qagiu (gsl_function * f,
  156.                double a,
  157.                double epsabs, double epsrel, size_t limit,
  158.                gsl_integration_workspace * workspace,
  159.                double *result, double *abserr)
  160. {
  161.   int status;
  162.  
  163.   gsl_function f_transform;
  164.   struct iu_params transform_params  ;
  165.  
  166.   transform_params.a = a ;
  167.   transform_params.f = f ;
  168.  
  169.   f_transform.function = &iu_transform;
  170.   f_transform.params = &transform_params;
  171.  
  172.   status = qags (&f_transform, 0.0, 1.0, 
  173.                  epsabs, epsrel, limit,
  174.                  workspace,
  175.                  result, abserr,
  176.                  &gsl_integration_qk15);
  177.  
  178.   return status;
  179. }
  180.  
  181. static double 
  182. iu_transform (double t, void *params)
  183. {
  184.   struct iu_params *p = (struct iu_params *) params;
  185.   double a = p->a;
  186.   gsl_function * f = p->f;
  187.   double x = a + (1 - t) / t;
  188.   double y = GSL_FN_EVAL (f, x);
  189.   return (y / t) / t;
  190. }
  191.  
  192. /* Main integration function */
  193.  
  194. static int
  195. qags (const gsl_function * f,
  196.       const double a, const double b,
  197.       const double epsabs, const double epsrel,
  198.       const size_t limit,
  199.       gsl_integration_workspace * workspace,
  200.       double *result, double *abserr,
  201.       gsl_integration_rule * q)
  202. {
  203.   double area, errsum;
  204.   double res_ext, err_ext;
  205.   double result0, abserr0, resabs0, resasc0;
  206.   double tolerance;
  207.  
  208.   double ertest = 0;
  209.   double error_over_large_intervals = 0;
  210.   double reseps = 0, abseps = 0, correc = 0;
  211.   size_t ktmin = 0;
  212.   int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
  213.   int error_type = 0, error_type2 = 0;
  214.  
  215.   size_t iteration = 0;
  216.  
  217.   int positive_integrand = 0;
  218.   int extrapolate = 0;
  219.   int disallow_extrapolation = 0;
  220.  
  221.   struct extrapolation_table table;
  222.  
  223.   /* Initialize results */
  224.  
  225.   initialise (workspace, a, b);
  226.  
  227.   *result = 0;
  228.   *abserr = 0;
  229.  
  230.   if (limit > workspace->limit)
  231.     {
  232.       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
  233.     }
  234.  
  235.   /* Test on accuracy */
  236.  
  237.   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
  238.     {
  239.       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
  240.          GSL_EBADTOL);
  241.     }
  242.  
  243.   /* Perform the first integration */
  244.  
  245.   q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
  246.  
  247.   set_initial_result (workspace, result0, abserr0);
  248.  
  249.   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
  250.  
  251.   if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
  252.     {
  253.       *result = result0;
  254.       *abserr = abserr0;
  255.  
  256.       GSL_ERROR ("cannot reach tolerance because of roundoff error"
  257.          "on first attempt", GSL_EROUND);
  258.     }
  259.   else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
  260.     {
  261.       *result = result0;
  262.       *abserr = abserr0;
  263.  
  264.       return GSL_SUCCESS;
  265.     }
  266.   else if (limit == 1)
  267.     {
  268.       *result = result0;
  269.       *abserr = abserr0;
  270.  
  271.       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
  272.     }
  273.  
  274.   /* Initialization */
  275.  
  276.   initialise_table (&table);
  277.   append_table (&table, result0);
  278.  
  279.   area = result0;
  280.   errsum = abserr0;
  281.  
  282.   res_ext = result0;
  283.   err_ext = GSL_DBL_MAX;
  284.  
  285.   positive_integrand = test_positivity (result0, resabs0);
  286.  
  287.   iteration = 1;
  288.  
  289.   do
  290.     {
  291.       size_t current_level;
  292.       double a1, b1, a2, b2;
  293.       double a_i, b_i, r_i, e_i;
  294.       double area1 = 0, area2 = 0, area12 = 0;
  295.       double error1 = 0, error2 = 0, error12 = 0;
  296.       double resasc1, resasc2;
  297.       double resabs1, resabs2;
  298.       double last_e_i;
  299.  
  300.       /* Bisect the subinterval with the largest error estimate */
  301.  
  302.       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
  303.  
  304.       current_level = workspace->level[workspace->i] + 1;
  305.  
  306.       a1 = a_i;
  307.       b1 = 0.5 * (a_i + b_i);
  308.       a2 = b1;
  309.       b2 = b_i;
  310.  
  311.       iteration++;
  312.  
  313.       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
  314.       q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
  315.  
  316.       area12 = area1 + area2;
  317.       error12 = error1 + error2;
  318.       last_e_i = e_i;
  319.  
  320.       /* Improve previous approximations to the integral and test for
  321.          accuracy.
  322.  
  323.          We write these expressions in the same way as the original
  324.          QUADPACK code so that the rounding errors are the same, which
  325.          makes testing easier. */
  326.  
  327.       errsum = errsum + error12 - e_i;
  328.       area = area + area12 - r_i;
  329.  
  330.       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
  331.  
  332.       if (resasc1 != error1 && resasc2 != error2)
  333.     {
  334.       double delta = r_i - area12;
  335.  
  336.       if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
  337.         {
  338.           if (!extrapolate)
  339.         {
  340.           roundoff_type1++;
  341.         }
  342.           else
  343.         {
  344.           roundoff_type2++;
  345.         }
  346.         }
  347.       if (iteration > 10 && error12 > e_i)
  348.         {
  349.           roundoff_type3++;
  350.         }
  351.     }
  352.  
  353.       /* Test for roundoff and eventually set error flag */
  354.  
  355.       if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
  356.     {
  357.       error_type = 2;    /* round off error */
  358.     }
  359.  
  360.       if (roundoff_type2 >= 5)
  361.     {
  362.       error_type2 = 1;
  363.     }
  364.  
  365.       /* set error flag in the case of bad integrand behaviour at
  366.          a point of the integration range */
  367.  
  368.       if (subinterval_too_small (a1, a2, b2))
  369.     {
  370.       error_type = 4;
  371.     }
  372.  
  373.       /* append the newly-created intervals to the list */
  374.  
  375.       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
  376.  
  377.       if (errsum <= tolerance)
  378.     {
  379.       goto compute_result;
  380.     }
  381.  
  382.       if (error_type)
  383.     {
  384.       break;
  385.     }
  386.  
  387.       if (iteration >= limit - 1)
  388.     {
  389.       error_type = 1;
  390.       break;
  391.     }
  392.  
  393.       if (iteration == 2)    /* set up variables on first iteration */
  394.     {
  395.       error_over_large_intervals = errsum;
  396.       ertest = tolerance;
  397.       append_table (&table, area);
  398.       continue;
  399.     }
  400.  
  401.       if (disallow_extrapolation)
  402.     {
  403.       continue;
  404.     }
  405.  
  406.       error_over_large_intervals += -last_e_i;
  407.  
  408.       if (current_level < workspace->maximum_level)
  409.     {
  410.       error_over_large_intervals += error12;
  411.     }
  412.  
  413.       if (!extrapolate)
  414.     {
  415.       /* test whether the interval to be bisected next is the
  416.          smallest interval. */
  417.  
  418.       if (large_interval (workspace))
  419.         continue;
  420.  
  421.       extrapolate = 1;
  422.       workspace->nrmax = 1;
  423.     }
  424.  
  425.       if (!error_type2 && error_over_large_intervals > ertest)
  426.     {
  427.       if (increase_nrmax (workspace))
  428.         continue;
  429.     }
  430.  
  431.       /* Perform extrapolation */
  432.  
  433.       append_table (&table, area);
  434.  
  435.       qelg (&table, &reseps, &abseps);
  436.  
  437.       ktmin++;
  438.  
  439.       if (ktmin > 5 && err_ext < 0.001 * errsum)
  440.     {
  441.       error_type = 5;
  442.     }
  443.  
  444.       if (abseps < err_ext)
  445.     {
  446.       ktmin = 0;
  447.       err_ext = abseps;
  448.       res_ext = reseps;
  449.       correc = error_over_large_intervals;
  450.       ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
  451.       if (err_ext <= ertest)
  452.         break;
  453.     }
  454.  
  455.       /* Prepare bisection of the smallest interval. */
  456.  
  457.       if (table.n == 1)
  458.     {
  459.       disallow_extrapolation = 1;
  460.     }
  461.  
  462.       if (error_type == 5)
  463.     {
  464.       break;
  465.     }
  466.  
  467.       /* work on interval with largest error */
  468.  
  469.       reset_nrmax (workspace);
  470.       extrapolate = 0;
  471.       error_over_large_intervals = errsum;
  472.  
  473.     }
  474.   while (iteration < limit);
  475.  
  476.   *result = res_ext;
  477.   *abserr = err_ext;
  478.  
  479.   if (err_ext == GSL_DBL_MAX)
  480.     goto compute_result;
  481.  
  482.   if (error_type || error_type2)
  483.     {
  484.       if (error_type2)
  485.     {
  486.       err_ext += correc;
  487.     }
  488.  
  489.       if (error_type == 0)
  490.     error_type = 3;
  491.  
  492.       if (res_ext != 0.0 && area != 0.0)
  493.     {
  494.       if (err_ext / fabs (res_ext) > errsum / fabs (area))
  495.         goto compute_result;
  496.     }
  497.       else if (err_ext > errsum)
  498.     {
  499.       goto compute_result;
  500.     }
  501.       else if (area == 0.0)
  502.     {
  503.       goto return_error;
  504.     }
  505.     }
  506.  
  507.   /*  Test on divergence. */
  508.  
  509.   {
  510.     double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
  511.  
  512.     if (!positive_integrand && max_area < 0.01 * resabs0)
  513.       goto return_error;
  514.   }
  515.  
  516.   {
  517.     double ratio = res_ext / area;
  518.  
  519.     if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
  520.       error_type = 6;
  521.   }
  522.  
  523.   goto return_error;
  524.  
  525. compute_result:
  526.  
  527.   *result = sum_results (workspace);
  528.   *abserr = errsum;
  529.  
  530. return_error:
  531.  
  532.   if (error_type > 2)
  533.     error_type--;
  534.  
  535.  
  536.  
  537.   if (error_type == 0) 
  538.     {
  539.       return GSL_SUCCESS;
  540.     }
  541.   else if (error_type == 1)
  542.     {
  543.       GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
  544.     }
  545.   else if (error_type == 2)
  546.     {
  547.       GSL_ERROR ("cannot reach tolerance because of roundoff error",
  548.          GSL_EROUND);
  549.     }
  550.   else if (error_type == 3)
  551.     {
  552.       GSL_ERROR ("bad integrand behavior found in the integration interval",
  553.          GSL_ESING);
  554.     }
  555.   else if (error_type == 4)
  556.     {
  557.       GSL_ERROR ("roundoff error detected in the extrapolation table",
  558.          GSL_EROUND);
  559.     }
  560.   else if (error_type == 5)
  561.     {
  562.       GSL_ERROR ("integral is divergent, or slowly convergent",
  563.          GSL_EDIVERGE);
  564.     }
  565.   else
  566.     {
  567.       GSL_ERROR ("could not integrate function", GSL_EFAILED);
  568.     }
  569.  
  570. }
  571.